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We study corrections to single tetrahedron based approximations for the entropy, specific heat and 
uniform susceptibility of the pyrochlore lattice Ising antiferromagnet, by a Numerical Linked Cluster 
(NLC) expansion. In a tetrahedron based NLC, the first order gives the Pauling residual entropy 
of | log | « 0.20273. A 16-th order NLC calculation changes the residual entropy to 0.205507 
a correction of 1.37 percent over the Pauling value. At high temperatures, the accuracy of the 
calculations is verified by a high temperature series expansion. We find the corrections to the single 
tetrahedron approximations to be at most a few percent for all the thermodynamic properties. 

PACS numbers: 74.70.-b,75.10.Jm,75.40.Gb,75.30.Ds 



INTRODUCTION 

The residual entropy of ice was first calculated in a 
classic paper by Pauling. [1] This entropy arises from the 
fact that an oxygen atom in ice is surrounded by four pro- 
tons. Any two of them can move in close to it, making 
up the H2O unit, while the other two stay away becom- 
ing part of the neighboring water molecules. Extensive 
ground state entropy is a hallmark of highly frustrated 
systems, and is well appreciated in context of Ising mod- 
els, at least since the exact solution of the triangular- 
lattice antiferromagnet. A direct connection between 
the entropy of ice and ground state entropy of Ising spin 
models was first shown by Anderson. [3[ Accurate calcula- 
tions of thermodynamic properties of systems with exten- 
sive ground state degeneracy [i|-[6[ remains a challenging 
task. 

Interest in such systems has grown considerably 
with the discovery of spin-ice materials. Q These are 
pyrochlore-lattice spin systems with strong uniaxial 
anisotropy. In each tetrahedron of the lattice, two 
spins point in and two point out, providing an ex- 
act realization of Anderson's spin-ice mapping. In the 
real materials, the spin-ice states arise from a combina- 
tion of nearest-neighbor exchange and long-range dipo- 
lar interactions. 0, Q Due to the unusual angles between 
the easy-axis directions at neighboring sites, antifcrro- 
magnetic exchange leads to a lower energy for all-in/all- 
out configurations in a tetrahedron, while ferromagnetic 
exchange leads to the two- in/two-out spin-ice configura- 
tions. The measured residual entropy, in these systems, 
is in good agreement with the Pauling value. Q 

In this paper, our focus is the study of an Ising antifer- 
romagnet on the pyrochlore lattice, which from a statis- 
tical mechanics point of view is equivalent to a ferromag- 
netically exchange coupled spin-ice material. It is defined 



by the Hamiltonian 



(1) 
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Here the spins Si take values ±1 and the sum runs over 
all nearest-neighbor bonds of the pyrochlore lattice. Sin- 
gle tetrahedron based approximations for the thermody- 
namic properties of this model are common, [H 13. [HI. [H, Ho. 
ll| We note that single tetrahedron based approximation 
for residual entropy is not the same as the entropy of a 
single tetrahedron. Translating Pauling's argument to 
the spin language, each spin has two states, but in each 
tetrahedron only 6 out of 16 spin configurations obey the 
ice-rules. Treating the constraint in each tetrahedron 
as independent, a system with N spins has 2 N (^) Nt 
ground states, where N t is number of tetrahedra. Since, 
N t equals N/2, this leads to a ground state entropy per 
spin of S = \ log | . 

Corrections to the Pauling expression for residual 
entropy of ice have been studied before. [l2| Follow- 
ing earlier work by Takahashi(l3| and DiMargio and 
Stilling er. I14l| Nagle used graphical methods on Vertex 
models 12] to study ground state residual entropy for cu- 
bic ice corresponding to the pyrochlore lattice as well as 
the more common ice structure known as Layered Hexag- 
onal ice. He found that the difference in residual entropy 
between the two structures was negligible. The correc- 
tions to the Pauling expression was estimated to be about 
1.1 percent. Comparable corrections (approximately 1.2 
percent) have been estimated in more recent computa- 
tional studies. .15] Here, we will calculate the thermody- 
namic properties of the antiferromagnetic Ising model on 
the pyrochlore lattice using series expansion methods. [3] 

In a Numerical Linked Cluster expansion (NLC), [13] 
an extensive property P for a large lattice C of N-sites is 
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expressed as 



P(C)/N = L(c) x W P {c) 



(2) 



The sum is over distinct clusters of the lattice. L(c) is 
the lattice constant of the cluster, given by the number 
of embeddings of the cluster in the lattice, per site. The 
quantity Wp(c) is the weight of the cluster associated 
with the property P, which is defined by the subgraph 
subtraction scheme 



W P (c)=P(c)-J2Wp(s), 



(3) 



where the sum runs over all subclusters of the cluster c. 
Thus, to carry out the calculation up to some order, one 
needs a count of all the clusters up to some order, and 
the property P(c) needs to be calculated for every cluster 
to high precision using numerical methods. 

The number of clusters, needed in the study, is signif- 
icantly reduced if it can be shown that only star-graphs 
contribute. 18[ Consider a graph as a collection of sites 
that are connected pair-wise by the bonds of the graph. 
A site is a point of articulation in a graph, if cutting 
all the bonds incident on the site, makes the rest of the 
graph disconnected. A star-graph is one that has no such 
articulation site. Star graph expansions have been used 
to develop high temperature expansions for various clas- 
sical spin models. [lSj But, they can also be used in a Nu- 
merical Linked Cluster (NLC) scheme, where no small 
expansion parameter is needed. (l7j Rather, thermody- 
namic properties of finite clusters can be calculated at 
a given temperature and then the principle of inclusion 
and exclusion can be used to calculate the thermody- 
namic property of the infinite system by summing up 
contributions from all allowed clusters. 

A star graph expansion requires that all articulated 
graphs have zero weight. This is very simple to see for 
properties that can be obtained from the logarithm of the 
zero-field partition function. For any articulated graph, 
one can show that a partition function, that is normalized 
to unity in the absence of interactions, becomes a product 
of partition functions over the two subgraphs articulated 
at a point. Hence \ogZ becomes a sum. Thus for an 
articulated graph c made up of parts a and b, which only 
share a point of articulation, the property satisfies 



P(c) = P(a) + P(b). 



(4) 



This is enough to ensure that after subgraph subtraction 
the articulated graph has zero weight. It is a little harder 
to see how to develop star-graph expansion for the uni- 
form susceptibility. One needs to consider the matrix M 
with elements (H,^ 
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FIG. 1: (color online) First two star-graphs of the pyrochlore 
lattice made up of complete tetrahedra. The first one is a 
single tetrahedron. The second is a ring of six tetrahedra 
that alternately point in and out of the board (denoted by + 
and — signs). 



where i and j are sites of the cluster and angular brackets 
denote thermal expectation values. One can show that 
for a graph articulated at site k into parts a and b (with 



spin normalization < St >= 1) 

< SiSj >=< SiSk > a < SkSj >b, 



(6) 



when i and j belong to a and b respectively. And within 
a and b, the correlations remain the same as in the sub- 
cluster. This can be used to further show that 19|, 2(1 
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(7) 



where Af _1 is the inverse of the matrix M, has a star 
graph expansion. Furthermore, for the infinite lattice, 
this is related to the inverse of the uniform structure 
factor, or temperature times the uniform susceptibility 
by the relation, 



i,3 



(8) 



Mi. 



(5) 



In a lattice of corner-sharing tetrahedra, it is natu- 
ral to consider an NLC expansion scheme, in which all 
interactions are grouped into tetrahedra. Thus, apart 
from a single site, all graphs are made up of complete 
tetrahedra. This scheme is particularly useful because 
the tetrahedra are joined at corners in the lattice. In a 
graph with two tetrahedra, the site where they are joined 
becomes a point of articulation. Cutting all the bonds 
at that site makes the graph disconnected. Thus, in a 
star-graph expansion, the two-tetrahedron graph makes 
no contribution. It also means that the single tetrahe- 
dron approximation becomes exact on a Husimi tree of 
tetrahedra. (21. 22 1 where there are no other closed loops 
of tetrahedra. 

The tetrahedra of a pyrochlore lattice are known to 
form a diamond lattice. Thus any graph counting prob- 
lem involving tetrahedra is equivalent to counting graphs 
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on the diamond lattice, where sites of the diamond lattice 
represent tetrahedra while bonds of the diamond lattice 
represent shared sites between neighboring tetrahedra. 
All star graphs on the diamond lattice up to 16 bonds 
have been listed by Sykes et aZ.[23j We will make use of 
these to calculate the expansions to 16-th order. One 
should note that all lattice constants of the diamond lat- 
tice need to be divided by 2, because the number of tetra- 
hedra is \ the number of sites in the pyrochlore lattice. 

We illustrate how the method works by showing the 
first two orders of the calculation for the ground state 
entropy (with the Boltzmann constant ks = !)• The 
first two star-graphs with complete tetrahedra are shown 
in Fig. 1. They have a lattice constant per site of the 
pyrochlore lattice of \ and 1 respectively. One also needs 
to consider a single-site, which provides all contributions 
before any interactions are included. It has a count of 
unity. 

For our illustration, the ground state entropy is the 
property P. In zeroth order, the single site has two 
ground states, giving an entropy of log 2. It has no sub- 
graphs. Hence, its weight is also log 2. The first star 
graph, a single tetrahedron, has 6 ground states. Hence 
the property P for the graph equals log 6. To obtain 
its weight, one must subtract the weights of the 4 sites. 
Thus the weight of the single tetrahedron is 

W = log6-41og2 = log3/8 

Thus, to first order the ground state entropy, per site, for 
the infinite system is 

S = log 2 + - log 3/8 = - log 3/2 w 0.20273. 

Note that the factor of | in front of log 3/2 is the count 
of number of tetrahedra per lattice site, which is one- 
half. This is the Pauling answer. The next star graph is 
a graph of 6 tetrahedra (See Fig. 1). It has 730 ground 
states. It has 6 single tetrahedron subgraphs and 18 sites. 
Thus, its weight is 

W = log 730 - 6 log 3/8 - 18 log 2 = log 730/729 
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FIG. 2: (color online) Entropy of the pyrochlore lattice Ising 
antiferromagnet as a function of temperature. 
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FIG. 3: (color online) Heat capacity of the pyrochlore lattice 
Ising antiferromagnet as a function of temperature. 



Thus, to this order, the entropy per site, becomes 



log2 + ilog3/8 



log 730/729 « 0.20410. 



These first corrections are analogous to triangle-based 
NLC calculations done by Rigol et al for the kagome 
lattice. ITJ In that case, the Pauling expression for 
ground state entropy is 0.50136. Adding the next cor- 
rection brought the entropy to 0.50182, much closer to 
known exact answer of approximately 0.50183. [li) Be- 
cause the count of longer loops increases much more 
rapidly in a 3-dimensional lattice, a high order calcula- 
tion is needed to assess the corrections more accurately. 



Here we have done only the first corrections for the uni- 
form susceptibility, but a 16-th order correction for the 
entropy and specific heat. 

High temperature series expansions for this system 
were derived some years ago: to order 19 for log Z (See 
Ref. 25 1) and to order 16 for X-[26| We have reanalyzed 
these series using Pade and d-log Pade approximants to 
yield the series results for comparison at high tempera- 
tures. 

Fig. 2 shows a plot of the entropy. The single tetra- 
hedra approximation is denoted NLC-1. NLC-6 includes 
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FIG. 4: (color online) Uniform susceptibility of the pyrochlore 
lattice Ising antiferromagnet as a function of temperature. 
Note that NLC-1 and NLC-6 can barely be differentiated in 
the plot. 



FIG. 5: (color online) Percentage correction to the quantities, 
entropy (P s ), heat capacity (P c ) and susceptibility (P x ) after 
6th and 16th order NLC calculations (denoted by the numbers 
6 and 16 in the legends). 



the next order correction. NLC-16 gives the result up to 
16th order. Fig. 3 shows the corresponding plots of the 
heat capacity. Fig. 4 shows plots of the uniform suscepti- 
bility. The high temperature expansions converge really 
well only above a temperature of 2. Note that the spins 
are normalized to be ±1. The temperature scale would 
be four times lower if they were normalized to ±1/2. 

In all cases, the single tetrahedron based approxima- 
tion is quite accurate. Let us define percent correction 
for the three quantities a — entropy (s), heat capacity 
(c) and susceptibility (x) as 

P a = 100 * {a{N) - a(l))/a(l). (9) 

Here a(N) is the quantity after iVth order NLC. and a(l) 
is the property with only the one tetrahedron cluster. A 
plot of P a with temperature is shown in Fig. 5. The 
largest corrections are near T — 0. They are a little over 
one percent for the entropy and about 5 percent for the 
specific heat. 

At T — 0, our graphical scheme, reduces to one of 
counting ground states on increasingly larger clusters. 
While this scheme is not identical to the one used by 
Nagle,[12j it is closely related. As found by Nagle, the 
first few orders give the same answer for the cubic and 
layered hexagonal structures of ice. We also find that the 
first corrections to the Pauling entropy are identical for 
the two lattices. Since, our interest is in the pyrochlore 
lattice spin model, we have not studied the layered hexag- 
onal structure in higher orders. 

In conclusion, we have used a star-graph expansion 
scheme to show that the Pauling approximation for the 



entropy of ice and calculation of other properties of the 
pyrochlore lattice Ising model based on a single tetrahe- 
dron is really a first term in a Numerical Linked Cluster 
(NLC) scheme. In the corner sharing lattice, this scheme 
is highly accurate. We have calculated corrections to the 
single tetrahedron approximations and found them to be 
at most a few percent for different thermodynamic quan- 
tities. 

These ideas of star-graph expansions are also valid 
for continuous spin models, where also single tetrahedra 
based schemes have been used successfully, [lfj However, 
it should be noted that finite-size calculations of proper- 
ties of the continuous spin systems can be a non-trivial 
task. Star graph expansions are also valid in presence of 
dilution and quenched disorder. [20ll27l | Star graph expan- 
sions are not valid for quantum spin models and hence 
all connected clusters of tetrahedra need to be included 
in these cases. 
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